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Abstract 

It has been known that the space-time CE/SE method can be used to obtain ID, 2D, and 3D steady and 
unsteady flow solutions with Mach numbers ranging from 0.0028 to 10. However, it is also known that a CE/SE 
solution may become overly dissipative when the Mach number is very small. As an initial attempt to remedy this 
weakness, new ID Courant number and Mach number insensitive CE/SE Euler solvers are developed using several 
key concepts underlying the recent successful development of Courant number insensitive CE/SE schemes. 
Numerical results indicate that the new solvers are capable of resolving crisply a contact discontinuity embedded in 
a flow with the maximum Mach number = 0.01. 


1. Intro dution 

The space-time conservation element and solution element (CE/SE) method is a high-resolution and 
genuinely multidimensional method for solving conservation laws [1-58]. Its nontraditional features include: 
(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution 
elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a novel time marching strategy 
that has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated 
without using any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability 
to capture shocks without using Riemann solvers); (iv) the requirement that each scheme be built from a 
non-dissipative core scheme and, as a result, the numerical dissipation can be controlled effectively; and (v) 
the fact that mesh values of the physical dependent variables and their spatial derivatives are considered as 
independent marching variables to be solve for simultaneously. Note that CEs are nonoverlapping space-time 
subdomains introduced such that (i) the computational domain can be filled by these subdomains; and (ii) 
flux conservation can be enforced over each of them and also over the union of any combination of them. 
On the other hand, SEs are space-time subdomains introduced such that (i) the boundary of each CE can 
be divided into several component parts with each of them belonging to a unique SE; and (ii) within a SE, 
any physical flux vector is approximated using simple smooth functions. In general, a CE does not coincide 
with a SE. 

Without using flux-splitting or other special techniques, since its inception [1] the unstructured-mesh 
compatible CE/SE method has been used to obtain numerous accurate ID, 2D and 3D steady and unsteady 
flow solutions with Mach numbers ranging from 0.0028 to 10 [42]. The phyical phenomena modeled include 
traveling and interacting shocks, acoustic waves, shedding vortices, viscous flows, detonation waves, cavi- 
tation, flows in fluid film bearings, heat conduction with melting and/or freezing, electrodynamics, MHD 
vortex, hydraulic jump, crystal growth, and chromatographic problems [2-58]. In particular, the rather 
unique capability of the CE/SE method to resolve both strong shocks and small disturbances (e.g., acoustic 
waves) simultaneously [11,13,14] makes it an effective tool for attacking computational aeroacoustics (CAA) 
problems. Note that the fact that second-order CE/SE schemes can solve CAA problems accurately is an 
exception to the commonly-held belief that a second-order scheme is not adequate for solving CAA problems. 
Also note that, while numerical dissipation is needed for shock capturing, it may also result in annihilation of 
small disturbances. Thus a solver that can handle both strong shocks and small disturbances simultaneously 
must be able to overcome this difficulty. 

In spite of its past successes, there is still room for improving the CE/SE method. An example is 
the fact that, in a CE/SE simulation with a fixed total marching time, generally the numerical dissipation 
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increases as the value of the Courant number decreases from 1, its maximum stability bound. As such, in a 
case with a large Courant number disparity (e.g., a simulation with a highly non-uniform spatial mesh and a 
spatially independent time step), the performance sensitivity with respect to the Courant number can lead 
to a solution that is highly dissipative in a region where the local Courant number <C 1. Another example 
is the fact that a CE/SE solution may become overly dissipative when Mach number is very small. 

By using the fact that each CE/SE scheme is built from a non-dissipative core scheme, robust one- 
dimensional and multidimensional Courant number insensitive schemes were described in [50,52]. As a follow- 
up, new one-dimensional Courant number and Mach number insensitive CE/SE Euler solvers have been 
developed using several key concepts underlying the development of Courant number insensitive schemes. 
Numerical results indicate that the new solvers are capable of resolving crisply a contact discontinuity 
embeded in a flow with the maxmum Mach number = 0.01. The rest of the paper is organized as follows. 
A review of the existing CE/SE schemes is given in Secs. 2-4. The new Courant number and Mach number 
insensitive schemes are described in Sec. 5. Numerical results are presented in Sec. 6. Conclusions and 
discussions are given in Sec. 7. 

2. Review of the basic ID CE/SE method 

For simplicity, we review the existing CE/SE schemes for the PDE 


du du 

m +a d^ = 0 


(2.1) 


where a ^ 0 is a constant. Let X\ = x, and X2 = t be considered as the coordinates of a two-dimensional Eu- 
clidean space E^. Then, because Eq. (2.1) can be expressed as V • h = 0 with h = (au, u), Gauss’ divergence 
theorem in the space-time E 2 implies that Eq. (2.1) is the differential form of the integral conservation law 


(f h-ds = 0 (2.2) 

Js{v ) 

As depicted in Fig. 1, here (i) S/E) is the boundary of an arbitrary space-time region V in E 2 , and (ii) 
ds = dan with da and ft, respectively, being the area and the unit outward normal of a surface element 
on S(V). Note that: (i) because h ■ ds is the space-time flux of h leaving the region V through the surface 
element ds, Eq. (2.2) simply states that the total space-time flux of h leaving V through S(V) vanishes; 
(ii) in E 2 , da is the length of a line segment on the simple closed curve S(E); and (iii) all mathematical 
operations can be carried out as though E 2 were an ordinary two-dimensional Euclidean space. 

To proceed, let ff denote the set of all space-time staggered mesh points in E 2 (dots in Fig. 2(a)), where 

n = 0, ±1/2, ±1, ±3/2, ±2, . . ., and, for each n, j = n± 1/2, n±3/2, n±5/2, Each (j, n) € is associated 

with a solution element, i.e., SE(j, n). By definition, SE (j, n) is the interior of the space-time region bounded 
by a dashed curve depicted in Fig. 2(b). It includes a horizontal line segment, a vertical line segment, and 
their immediate neighborhood. 

Let (x,t) € SE (j,n). Then Eq. (2.2) will be simulated numerically assuming that u(x,t) and h(x,t), 
respectively, are approximated by 


and 


u*(x,t;j,n ) = Uj + {u x )j(x - Xj) + (u t )"(f - t n ) (2.3) 

h*(x,t\j,n) = f ( au*(x,t\j,n ), u*(x,t\j,n )) (2.4) 


Note that (i) u”, (u x )", and (w t )" are constants in SE(j, n), (ii) ( Xj,t n ) are the coordinates of the mesh point 
(j, n ) with Xj = jAa: and t n = nAt, and (iii) Eq. (2.4) is the numerical analogue of the definition h = (au, u). 

Let u = u*(x,t\j,ri) satisfy Eq. (2.1) within SE(j, n). Then one has (u t )" = — a(u x )J. As a result, 
Eq. (2.3) reduces to 

u*(x,t;j,n ) = + (u x )j [(x - Xj) -a(t- t n )\ (2.5) 
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i.e., u'j and (u x )” are the only independent marching variables associated with ( j,n ). 

Let E ‘2 be divided into nonoverlapping rectangular regions (see Fig. 2(a)). As depicted in Figs. 2(c)-2(e), 
(i) two such regions, i.e., CE _(/, n) and CE+(j, n), are associated with each interior mesh point (j,n) £ Q: 
and (ii) CE(j, n) is the union of CE_(j, n) and CE + (j, n). 

Given the above preliminaries, we are ready to describe the existing CE/SE solvers for Eq. (2.1). 

2.1. The a scheme 

Note that, among the line segments forming the boundary of CE_(j, n), AB and AD belong to SE(j. n), 
while CB and CD belong to SE(j — 1/2, n — 1/2). Similarly, the boundary of CE + (j, n) belongs to either 
SE(j, n) or SE(j + 1/2, n — 1/2). As a result, by imposing two conservation conditions at each (j, n) € f 2, 

i.e., 

(f h*-ds = 0 and (f h* ■ ds = 0, (j,n) £ O (2.6) 

JS(CE+(j,n)) J S{CE- (j,n)) 

and using Eqs. (2.4) and (2.5), one can obtain two equations for the two unknowns it" and (u x )". In fact, 
let (i) v = f aAt/Ax , and (ii) for any (j,n) £ fl, 

(**)]=!{-(«*)] ( 2 - 7 ) 

then Eq. (2.6) implies that (i) 

U j = j{( 1 + V )Vl/2 + f 1 - ^) U "+l/2 2 + (! - ^) ~ ( U s)"+l/2 } ( 2 ' 8 ) 

and, assuming \u\ ^ 1, (ii) 

(«*)" = («!)" (M^l) (2-9) 

with 

04)" = f \ '«"+ 1/2 “ u j-: 1/2 - (! + ^ («*)"+] 1/2 - (! - v ){ u xTjZl/l (W\ ± 1) (2-10) 

The a scheme, i.e., the inviscicl case of the a-fj, scheme [1,3,9], is formed by Eqs. (2.8) and (2.9). Note that, 
because 

du ax du 

dx 4 dx 

if x = f x/(ax/4), the normalized parameter («$)" may be interpreted as the value at (j,n) of the derivative 
of u with respective to the normalized coordinate x. Also note that the superscript symbol “a” in the 
parameter (u|)” is introduced to remind the reader that Eq. (2.9) is valid for the a scheme. 

The review of the a scheme is concluded with the following remarks: 

(a) Even though it is introduced to model a single PDE (i.e., Eq. (2.1)) with a single dependent variable 
u, the a scheme is formed by two coupled discrete equations (i.e., Eqs. (2.8) and (2.9)) involving two 
independent numerical variables and (tt x )”. It is shown in [1] that Eqs. (2.8) and (2.9) are consistent 
with a pair of PDEs with one of them being Eq. (2.1). 

(b) The a scheme has the simplest stencil, i.e., a triangle with a vertex at the upper time level and the other 
two vertices at the lower time level. Furthermore, the number of the independent marching variables 
associated with a mesh point (j. n) £ ft is equal to the number of the mesh points at the (n — l/2)th 
time level that are part of the stencil. Note that the same relation also holds for many 2D and 3D 
CE/SE schemes [7,8,41]. 

(c) As shown in [3], the two amplification factors of the a scheme are identical to those of the leapfrog 
scheme. As a result, the a scheme is non-dissipative and it is stable if |^| < 1 (see the additional 
discussions given in Sec. 2.2). 

(d) Note that derivation of Eqs. (2.8) and (2.9) can be facilitated by the following observations: because 
u*(x,t'i j , n) is linear in x and t, it can be shown that the total flux of h* leaving CE_ (j, n) or CE + (j, n) 
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through any of the four line segments that form its boundary is equal to the scalar product of the vector 
h* evaluated at the midpoint of the line segment and the “surface” vector (i.e., the unit outward normal 
multiplied by the length) of the line segment. 

(e) Because, for any (j, n) G £2, the total flux of h* leaving each of CE_(j, n) and CE +(j,n) vanishes (see 
Eq. (2.6)), CE-(j,n ) and CE + (j,n), ( j,n ) € £2, will be referred to as the conservation elements (CEs) 
of the a scheme. In addition, because (i) the vector h* at any surface element lying on any interface 
separating two neighboring CEs is evaluated using the information from a single SE, and (ii) the unit 
outward normal vector on the surface element pointing outward from one of these two neighboring CEs 
is exactly the negative of that pointing outward from another CE, one concludes that the flux leaving one 
of these CEs through the interface is the negative of that leaving another CE through the same interface. 
As a result, the local conservation relations Eq. (2.6) lead to a global flux conservation relation, i.e., the 
total flux of h* leaving the boundary of any space-time region that is the union of any combination of 
CEs will also vanish. In particular, because CE(j, n) is the union of CE ~(j,n) and CE + (j, n), 

< j) h* ■ ds = 0, (j,n) £ £2 (2-11) 

JS(CE(j,n)) 

must follow from Eq. (2.6). In fact, it can be shown that Eq. (2.11) is equivalent to Eq. (2.8). 

(f) In addition to the non-dissipative a scheme, as will be shown, there is a family of its dissipative extensions 
in which only the less stringent conservation condition Eq. (2.11) is assumed [3]. Because Eq. (2.11) is 
equivalent to Eq. (2.8), for each of these extensions, n” is still evaluated using Eq. (2.8) while («$)" is 
evaluated using an equation different from Eq. (2.9). 

2.2. The a-e scheme and the c scheme 

To proceed, consider any (j, n) G £2. Then ( j ± 1/2, n — 1/2) G 12. Let 

«i±i/2 = <±1/2 + (**/2) (2-12) 


With the aid of Eq. (2.7) and the fact that the Courant number v = f aAt/Ax, a substitution of the relation 
(itt)" = — a(u x )(j into Eq. (2.12) results in 

<± 1/2 = (u-2vu S; ) n j ~l / / l (2.13) 

Note that, to simplify notation, in the above and hereafter we adopt a convention that can be explained 
using the expression on the right side of Eq. (2.13) as an example, i.e., 

(u - 2 - 2v(u S: ) n j ~l / / l 


Also note that, by definition, (j ± 1/2, n) ^ 12 if (j, n) G £2. Thus t^±i /2 is associated with a mesh point 
</ £2. The reader is warned that similar situations may occur in the rest of this paper. 

According to Eq. (2.12), can be interpreted as a first-order Taylor’s approximation of u at 

(j ± 1/2, n). Thus 


(«£)? = 


def AX / U j+l/2 U j—l/2 


AX 


(2.14) 


is a central-difference approximation of du/dx at (j. n). Note that: (i) the superscript “c” 
the reader of the central-difference nature of the term (u|)”; and (ii) by using Eqs. (2.13), 
one has 


KY, 



2 vux) 


ra-l/2 
3+ 1/2 


{u 


2 vuY) 


71-1/2' 

1 - 1/2 


is used to remind 
(2.14) and (2.10), 

(2.15) 
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and 


(4 


«£)? 



n- 1/2 
1 + 1/2 


+ ( u $) 


n- 1/2' 

J — 1/2 


The a-e scheme is formed by Eq. (2.8) and 


I iV 1 - 1 / 2 - 

4 { U j+ 1/2 M i-l/2 / 


(2.16) 


( Us )y = («?)7 + 2 C («|- U |)5 


(2.17) 


where e is a real number. Obviously the a-e scheme reduces to the a scheme when e = 0. Also, for the case 
e = 1/2, Eq. (2.17) reduces to 

(«*)? = («£)? (2.18) 

Because (it|)!/ represents a central-difference approximation, hereafter, to simplify its frequent references, 
the special a-e scheme with e = 1/2 will be referred to as the c scheme. 

To proceed, several key remarks about the a-e scheme are presented: 

(a) At each mesh point (j,ri) € Q, Eqs. (2.8) and (2.9) are the results of Eq. (2.6). Because Eq. (2.17) does 
not reduce to Eq. (2.9) except in the special case e = 0, at each mesh point ( j,n ) € f 1, generally the 
a-e scheme satisfies only the single conservation condition Eq. (2.11) (which is equivalent to Eq. (2.8)) 
rather than the two consevation conditions Eq. (2.6). However, because (u|)” generally is present on 
the right side of Eq. (2.17), the a-e scheme generally will still be burdened with the cost of solving two 
conservation conditions at each mesh point. The exception occurs only for the special case e = 1/2 
(i.e., the c scheme) in which Eq. (2.17) reduces to Eq. (2.18). As it turns out, implementation of 
a multidimensional Euler version of the c scheme does not require inverting any system of equations 
while a similar implementation involving a version of any other a-e scheme (e ^ 1/2) generally requires 
inverting, per mesh point and per time step, a system of several linear equations (to be exact, a system 
of eight and fifteen equations, respectively, for 2D and 3D Euler equations) [7,8,41]. As such, it is much 
more cost effective to use a multidimensional Euler version of the c scheme than using that of any other 
a-e scheme. Partly for this reason, extensions of the c scheme have been used extensively. 

(b) For the a-e scheme, it is shown in [3] that the principal and spurious amplification factors per At, 
respectively, are (A + ) 2 and (A_) 2 with 

A±(e, v , 9) = f ecos(0/2) — ivsm(6/2) ± \J (1 — e) [(1 — e)cos 2 (0/2) + (1 — i/ 2 )sin 2 (0/2)] (2.19) 

Here (i) i = f yf—l, and ( ii) 9, — 7 r < 9 < n, is the phase angle variation per ax. In addition, it is shown 
that (i) the necessary and sufficient conditions for the stability of the a-e scheme are 

0 < e < 1, and \u\ < 1 (2.20) 

and (ii) the a-e scheme becomes progressively diffusive as the value of e increases from 0 to 1. Note 
that, unless specified otherwise, in the remainder of the paper the ranges of e, v and 9, respectively, are 
defined by Eq. (2.20) and —tt < 9 < tt. 

(c) Let k be a constant. Then u = e lk ^ x ~ at ^ represents a plane wave solution to Eq. (2.1). For this solution 


the exact amplification factor per At 


def e ifc [ x-a ( t+At )l 
gik(x-at) 


g—ikaAt g— iv6 


( 2 . 21 ) 


where 9 = kAx. 

(d) According to Eq. (2.19), [A±(0, v, 9)] 2 , the amplification factors of the a scheme (which corresponds to 
the case e = 0) per At, have the following properties: 


|[A±(0, v, 0)] 2 | = 1 


( 2 . 22 ) 
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(2.23) 


iim [A + (°,^,6»)] 2 = e TlB 

Z /— >±1 

lim [A_(0,^,6»)] 2 = e ±id (2.24) 

1/— >±1 

and 

[A±(0, 0, 9)] 2 = 1 (2.25) 

On the other hand, e ~ w6 , the exact amplification factor per At, has the following properties: 

\e~ ive \ = l (2.26) 

lim e - ivf> = e^ iS (2.27) 

^—►±1 

and 

e~ ive = 1 if v = 0 (2.28) 

For the a scheme, Eqs. (2.22)-(2.28) imply that: (i) the two amplification factor of the scheme, and 
the exact amplification factor all have the same constant absolute value (= 1) and, thus, the scheme is 
non-dissipative; (ii) in the limit of |^| — » 1 (i.e. , v — > 1 or v — » —1), the principal amplification factor is 
identical to the exact amplification factor and, thus, the former has no dissipative or dispersive error in 
this limit; (iii) also in the limit of \v\ — ■> 1, the phase angle associated with the spurious amplification 
factor is exactly the negative of that associated with the exact amplification factor and, thus, the 
spurious amplification factor has a large dispersive error in this limit except when \9\ <C l(i.e., when 
the wavelengths of the errors 1); and (iv) when v = 0, the two amplification factors of the scheme, 
and the exact amplification factor are all equal to 1 and, thus, the two amplification factors of the 
scheme have no dissipative or dispersive error if v = 0. Because the accuracy of a scheme is essentially 
hinged on the behaviors of the principal amplification factor [1], according to the facts stated above, 
the a scheme tends to become very accurate when \v\ approaches 1 or 0. However, the short- wavelength 
errors associated with the spurious amplification factor (which could be introduced at t = 0 as a result of 
an inaccurate initial-value specification [1]) may appear in a solution as persistent (i.e., non-dissipative) 
numerical wiggles when \v\ approaches 1 [1,9]. 

(e) According to Eq. (2.19), [A± (1/2, v, 9)] 2 , the amplification factors of the c scheme (which corresponds 
to the case e = 1/2) per At, have the following properties: 


hm [A+(l/2,^, 6>)] 2 = e Tie 

Z '— »±1 


(2.29) 


lim [A_ (1/2, u, 9)] 2 = — sin 2 (0/2) 


(2.30) 


and 


[A±(l/2, 0, 9)] 2 


1 ' 
2 . 


1 ± cos (0/2) 1 / 2 " 


cos 2 (9/2) 


(2.31) 


For the c scheme, Eqs. (2.27)-(2.31) imply that: (i) in the limit of |i^| — > 1, the principal amplification 
factor is identical to the exact amplification factor and, thus, the former has no dissipative or dispersive 
errors in this limit; (ii) also in the limit of \v\ — > 1, the spurious amplification factor has large dissipative 
and dispersive errors; and (iii) when v = 0, both the principal and spurious amplification factors 
generally have large dissipative errors but no dispersive errors. According to the facts stated above, like 
the a scheme, the c scheme also tends to become very accurate when |^| approaches 1. However, unlike 
the a scheme, the errors associated with the spurious amplification factor of the c scheme generally do 
die out rapidly when |^| approaches 1. Also, in sharp contrast to the a scheme, the c scheme becomes 
highly dissipative when v approaches 0. 

From the above discussions, one concludes that: 
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(a) The advantages of the a scheme include: (i) it is non-dissipative throughout the range of e and v defined 
in Eq. (2.20); and (ii) when the value of \v\ is close to 0 or 1, the scheme is very accurate. On the other 
hand, its disadvantages include: (i) because it is non-dissipative, its extensions for nonlinear equations 
generally are unstable; (ii) when the value of \u\ is close to 1, the short- wavelength errors associated 
with the spurious amplification factors will not die out rapidly and, therefore, appear in a solution as 
persistent numerical wiggles; and (iii) comparing with the c scheme, it costs more to implement. 

(b) The advantages of the c scheme include: (i) when the value of \v\ is close to 1, it is very accurate and 
the short-wavelength errors associated with the spurious amplification factor also die out rapidly; (ii) 
because it is dissipative, its extensions for nonlinear equations generally are stable; and (iii) in terms of 
ease of implementation and computer cost, it is much more superior than any other a-e scheme. On the 
other hand, the c scheme has a serious disadvantage, i.e., it is very dissipative when v approaches 0. 


In Sec. 3, it will be shown that new solvers of Eq. (2.1) can indeed be constructed such that they possess all 
the advantages but none of the disadvantages listed above. Specifically, each of these solvers will be formed 
by Eq. (2.8) and a new equation in which (u $ )" is eva luted using a simple central-differencing procedure 
similar to that used to obtain (u|)”. In addition, (w $ )™ so obtained will be (i) identical to (u§)" in the limit 
of \u\ — > 1, and (ii) identical to (u|)" when v = 0. As such, each of these new solvers (i) is comparable to 
the c scheme in ease of implementation; (ii) becomes the c scheme in the limit of \v\ — * 1; and (iii) becomes 
the a scheme when v = 0. 


2.3. The w-a scheme — a special wiggle-suppressing scheme 


If discontinuities are present in a numerical solution, any a-e scheme such as the c scheme is not equipped 
to suppress numerical wiggles that generally appear near these discontinuities. To serve as a preliminary 
for future development, here we shall briefly review an extension of the c scheme which was introduced as a 
remedy for this deficiency [3,35]. 

To proceed, let 


(«*-)? = -r 


def AX [ u j u j- 1/2 


ax/2 


(2.32) 


and 


(«*+)? = 


def AX / u j+ 1/2 ~ u j 


ax/2 


(2.33) 


i.e., (ux-)(j and (w $+ )" are numerical analogues of du/dx at (j,n) 
respectively. It can be shown that 


(<4)j = \ (Ux- + %+)" 


evaluated from the left and the right, 


(2.34) 


i.e., (u|)" is the simple average of ('«$-)" and (w s+ )". As such, the c scheme can be extended by replacing 
(u§)" in Eq. (2.18) with an weighted average of (?%,_)" and (u s+ )". In other words, the resulting extension 
is formed by Eq. (2.8) and 

(%)” = («>-)?(«*-)? + (w+)](u s+ )] (2.35) 

where («;_)" and ( w+)' ), the weight factors associated with ( u s _)" and (u$+)" respectively, must satisfy the 
condition 

(«,_)” + («;+)? = 1 (2.36) 

at all (j, n) G O. In addition, the expression on the right side of Eq. (2.35) represents an interpolation (rather 
than an extrapolation) of (it^-)” and (u^+)j if and only if 

(w-)] > 0 and ( w + )? > 0 (2.37) 


For real variables x_, x+, and a > 0, let W- and W + be the functions defined by: (i) W-(x~, x + ; a) = 
W+(x_, x+; a ) = 1/2 if x_ = x+ = 0; and (ii) 


VF_(x_,x+; a) 


i^+r 

|x_|“ + |x+| a 


and W+ (x - , x+ ; a) 


\x-\ a + |x+| a 


(2.38) 
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(2.39) 


if either 7^ 0 or x+ 7^ 0. Furthermore, let 

(«; ± )7 = W t ((« iE _)7, («*+)?, a) 

Then (w_)" and (w + )” so defined satisfy Eqs. (2.36) and (2.37) and have the property that 

(«>-)" = («,+)? = 1/2 if a = 0 or | («$-)” | = | («*+)" I (2-40) 

Note that: (i) to avoid dividing by zero, in practice a small positive number such as 10~ 20 is added to each 
of the denominators in Eq. (2.38); and (ii) the special cases of Eq. (2.38) with a = 1 and a = 2 are used in 
the slope-limiter proposed by van Leer [59] , and van Albada et a 1. [60] . 

An extension of the c scheme is formed by Eqs. (2.8) and (2.35) with (tu_)” and (tt?+)” being defined 
by Eq. (2.39). Because it involves an weighted average which is dependent on a parameter a, hereafter the 
scheme is referred to as the w-a scheme. Let a > 0 and |(u s _)”| 7^ |(u s+ )"|. Then Eq. (2.38) implies that, 
of (ttjg_)” and (tt$+)", the one with smaller absolute value is associated with an weight factor >1/2. This 
observation coupled with Eqs. (2.34)-(2.37) leads to the conclusion that, of (uj_)" and (u i+ )", (w s )" will 
have an algebraic value closer to the one with smaller absolute value if (it^)™ is evaluated as an weighted 
average of (u S -)" and (% + )" according to Eq. (2.35). As a result, (u^)" so evaluated has a smaller absolute 
value than that evaluated using Eq. (2.18). In turn, numerical wiggles or overshoots can be annihilated by 
the additional numerical dissipation introduced as a result of this local “flattening” of (tt $ )™. It has been 
shown numerically that the extension is stable if \v\ < 1 and a > 0. Moreover, as a result of Eqs. (2.18), 
(2.34), (2.35) and (2.40), (i) the extension reduces to the c scheme when a = 0; and (ii) even if a > 0, 
the extension behaves very much like the c scheme in any smooth solution region (where the condition 
(its-)" = (ux+)'j more or less prevails) or at a solution extremum (where the condition (« s _)" = — (w$+)" 
more or less prevails). As such, the wiggle-suppressing power of the extension takes effect only if a > 0 and 
only in a solution region where |(u $ _)"| and |(ux+)”| differ substantially. 

3. The c-t and c-t* schemes 

In this section, the ideal solvers of Eq. (2.1) mentioned at the end of Sec. 2.2 will be constructed. As a 
preliminary, we shall show that («§)” can also be cast into a central-difference form when v = 0. 

To proceed, note that by assumption a/0. Thus v = 0 if and only if At = 0. Because \EF\ = |AD| = 
\CB\ = 0 (see Figs. 2(c,d)) when At = 0, the two conservation conditions given in Eq. (2.6) for the case 
v = 0, respectively reduce to the following conditions: (i) the flux leaving CE + (j,ri) through the top face 
AF is equal to that entering the same CE through the bottom face ED\ and (ii) the flux leaving CE_(j, n) 

through the top face AB is equal to that entering the same CE through the bottom face CD. According 

to Remark (d) given at the end of Sec. 2.1, the flux leaving CE + (j, n) through the top face AF is equal to 
the value of u* at the midpoint of AF (evaluated using the marching variables at point A) multiplied by 
|AF|, while that entering it through the bottom face ED is equal to the value of u* at the midpoint of ED 
(evaluated using the marching variables at point E) multiplied by \ED\. With the aid of these observations 
and the fact that \AF\ = \ED\, the above condition (i) implies that, when v = 0, the value of u* at the 
midpoint of AF evaluated using the marching variables at point A is equal to that at the midpoint of ED 
evaluated using the marching variables at point E. As such, the first conservation condition in Eq. (2.6) is 
equivalent to 

(u + 1 i s )j = (u - UxYj+ljl ( v = °) (3-1) 

if v = 0. Similarly, by using the above condition (ii), it can be shown that the second conservation condition 
in Eq. (2.6) is equivalent to 

(u - Uj)" = (u + Ux)™~s/2 ( u = °) (3-2) 

if v = 0. Because Eqs. (2.8) and (2.9) (which form the a scheme) are equivalent to Eq. (2.6) if \v\ 7^ 1, they 
must be equivalent to Eqs. (3.1) and (3.2) when v = 0. In fact, by substracting Eq. (3.2) from Eq. (3.1), one 
obtains Eq. (2.9) where (u|)” is the reduced form of Eq. (2.10) for the case v = 0, i.e., 

(«£)? = \ [(« - - (u + « $ );: i 1 / 2 2 ] {v = 0) (3.3) 
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Moreover, by summing over Eqs. (3.1) and (3.2), one has the reduced form of Eq. (2.8) for the case v = 0. 

With the aid of Eq. (3.3) and the facts that: (i) (u — and (u + UxYjZi /2 > respectively, represent 
an approximation of u at the midpoint of ED and that at the midpoint of CD (see Fig. 2(c,d)); and (ii) 
the distance between the two midpoints referred to above is ax/2, it becomes obvious that, for the special 
case v = 0, (it§)™ is indeed a central-difference approximation of du/dx at ( j,n — 1/2) (which coincides with 
(j, n) when v = 0). QED 

According to the above discussions, construction of the ideal solvers defined at the end of Sec. 2.2 is 
hinged on finding central-difference approximations for (it$)” such that each approximation (i) becomes 
(u|)” in the limit of \u\ — > 1, and (ii) reduces to the expression on the right side of Eq. (3.3) when v = 0. 
As a result of these observations, these new solvers can easily be constructed as the subschemes of the c-r 
scheme, a new class of CE/SE solvers for Eq. (2.1) to be described immediately. 

3.1. The c-r scheme 

To proceed, refer to Fig. 3. Here M + and M~ denote the midpoints of AF and AB, respectively. Also 
P + and P~ are two points on BF that satisfy the following conditions: (i) P + coincides with M + if and 
only if P~ coincides with M~; (ii) P + is to the right (left) of M + if and only if P~ is to the left (right) 
of A / - ; and (iii) |M+P+| = \M~P~\, i.e., M+P+ and M~P~ have the same length. In addition, let the 
parameter r be defined by: (i) r = 0 if P + coincides with Af + ; (ii) tax /A = |Af+P+| if P + is to the 
right of Af + ; and (iii) t ax /A = —\M + P + \ if P + is to the left of M + . Obviously, it follows from the above 
definitions that (i) r = 0 if P~ coincides with A f - ; (ii) r ax/A = \M~P~\ if P~ is to the left of A and 
(iii) r ax/A = — \M~P~\ if P~ is to the right of A 1~ . 

Moreover, let 

u(P + ) = f [u + (At/2)u t - (1 - / r)it 5 ;]”“ 1 1 / / 2 2 (3.4) 

and 

u(P~) = f [u + (At/2)u t + (1 - T)Ux]™Zl /2 (3-5) 

With the aid of Eq. (2.7), it is seen that u(P + ) is a first-order Taylor’s approximation of u at P + evaluated 
using the marching variables at point E, while u(P~) is a first-order Taylor’s approximation of u at P~ 
evaluated using the marching variables at point C. Also note that, by using the relation (itj)" = — a(it x )", 
Eqs. (3.4) and (3.5), respectively, can be simplified as 

u{P + ) = [u-(2u+l- r)ux]"+i/ 2 2 ( 3 - 6 ) 

and 

u(P ~ ) = [u ~(2l/-l + t)Ux]™Zi /2 ( 3 - 7 ) 

At this juncture, note that P + and P~ generally lie outside of SE(j+l/2, n— 1/2) and SE(j— 1/2, n— 1/2), 
respectively. Yet here, by definition, u(P + ) and u(P~) are evaluated as though P + and P~ lie within 
(j -F 1/2, n — 1/2) and (j — 1/2, n — 1/2) , respectively. At first glance, the current practice is inconsistent with 
a previously established rule. However, as explained by the reasons given below, the definition of u{P + ) and 
u(P~) is perfectly legitimate: 

(a) Recall that solution elements were introduced such that the boundary of a CE can be divided into 
several component parts with each of them belonging to a unique solution element. As such, the flux 
over a component part that belongs to a special solution element, say SE(j, n), can be unambiguously 
determined in terms of the marching variables at the mesh point ( j,n ). In other words, in related to 
evaluation of any flux conservation condition over any CE, Eqs. (2.3)-(2.5) can be applied only to a 
point (x,t) £ SE (j,n). 

(b) On the other hand, u(P + ) and u(P~) introduced here have nothing to do with flux evaluation. In fact, 
they will be used only in the construction of some numerical analogues of du/dx at ( j,n ). 
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To proceed, note that: (i) the mesh point (j.n) (i.e., point A depicted in Fie. 3) is the midpoint of 
P-P+, and (ii) \P-P+\ = (1 + r) ax/2. Thus 


def AX f u{P+)-u{P ) \ 
J 4 \ (1 + r) ax/2 ) 


(r^-1) 


(3.8) 


represents a central-difference approximation of du/dx at the mesh point (j, n). Thus the new scheme formed 
by Eq. (2.8) and 

(«x)? = (fix)? (3-9) 

represents a solver for Eq. (2.1). Because (i) (u s )" represents a central-difference approximation of (u s )", 
and (ii) the approximation is associated with the parameter r, hereafter the new scheme will be referred to 
as the c-t scheme. 

To explore the c-t scheme, note that Eqs. (3.6)-(3.8) can be combined to yield 


(fix)? = 


2(1 + r) 


{i 


U-(2v+l- t)Ux]™ + i/2 - \ u - ( 2u - 1 + r ) U x]?-l/2 } ( T ^ - 1 ) 


(3.10) 


Moreover, by using Eqs. (2.10), (2.16) and (3.10), one has 


(fix)? = («i)5 


2r 

1 + r 


(«* 


«£)? 


K 1 - T ) 

2(1 + r) 



n - 1/2 
1 + 1/2 


(« i ) 


n — 1/2 
1 - 1/2 


(3.11) 


By comparing Eq. (3.11) with (2.17), one concludes that the c-t scheme generally is different from the a-e 
scheme. In fact, a special case of the c-t scheme can be turned into that of the a-e scheme and vice versa if 
and only if either (i) r = 1 or (ii) v = 0. For the case r = 1, Eq. (3.11) implies that (u$)” = (u§)". In other 
words, the c scheme is the special case of the c-t scheme with r = 1, a fact that can also be deduced from 
the observation that the points P + and P~ depicted in Fig. 3, respectively, coincide with points F and B 
(i.e., the mesh points (j + 1/2, n) and (j — 1/2, n )) if r = 1. On the other hand, it is seen that, when v = 0, 
the c-r scheme become the a-e scheme with e = 2t/( 1 + r). In fact one can further deduce that c-t scheme 
reduces to the a scheme if and only if v = r = 0. 

Because the c-r scheme is formed by two rather complicated equations (i.e., Eqs. (2.8) and (3.9)) 
involving two parameters v and r, it were not expected that its von Neumann stability conditions could 
be cast into an explicit analytical form. However, it is shown rigorously in [55] that the c-t scheme is von 
Neumann stable if and only if 


v 2 < 1, t> t 0 (is 2 ), and (n 2 , r ) ± (1,1) 


(3.12) 


where 


Note that: 

(a) It can be shown that 
(hi) 


ro 


/ \ def i 
+o(s) = < 


if s = 0 


4- s-2y/2(2- s-s 2 ) 


:-l + Vl-2s + 5s 2 
2s 


if 0 < s < yj- 


if < 5 < 1 


(3.13) 


[55]: (i) t 0 (s ) is continuous at s = 0; (ii) r 0 (s) is consistently defined at s = 3/11; 


lim t'(s) = lim t'(s) = 121/90 

.. 3 - 3 + 


where t' 0 (s ) = dT 0 (s)/ds\ (iv) r 0 (s) is strictly montonically increasing in the interval 0 < s < 1; (v) 
r 0 (l) = 1; and (vi) 

s < t 0 (s) < y/s, 0 < s < 1 (3-14) 
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(b) For any given fixed value of \v\ < 1, the c-r scheme tends to become more dissipative as the value of r 
increases from t q (v 2 ), its minimum stable value. 

(c) The function r Q defined here is different from the function r 0 defined in [52]. In fact, for any v £ [0, 1], 
the value of t 0 (v 2 ) defined here is identical to that of r 0 (|n|) defined in [52]. Because the analytical 
n-T 0 (V 2 ) relation depicted in Fig. 4 here is virtually identical to the numerical ^-r(|t'|) relation depicted 
in Fig. 4 of [52], the analytical stability conditions given here are in complete agreement with those 
generated numerically and stated in Eq. (3.12) of [52]. 

With the above preliminaries, the ideal solvers of Eq. (2.1) defined at the end of Sec. 2.2 will be 
constructed in Sec. 3.2. 

3.2. The c-r* schemes 

The value of r used in the c-r scheme generally can be chosen independent of v. Here we will introduce 
a subset of the c-r scheme in which r is a function of v 2 for each member of this subset. As a preliminary, 
note that, by using the properties of r D (s) presented earlier, it can be shown that there exist infinitely many 
choices of a strictly monotonically increasing smooth function h(s), 0 < s < 1, which possesses the following 
properties: 

/i(0) = 0; lim h(s) = 1; and h(s) > r 0 (s) if 0 < s < 1 (3.15) 

S — > 1 — 

For each function h(s ) that satisfies the above conditions, each member of the subset referred to earlier is 
defined using the relation 

r = h{v 2 ) ( v 2 < 1) (3.16) 

Note that, using the definition of h and Eq. (3.16), one can easily infer from Fig. 3 a simple relation between 
the value of v 2 and the locations of P + and P~ , i.e., as the value of v 2 increases from 0 to 1, P + will move 
away from M + and edge toward the mesh point (j + 1/2, n) while P~ will move away from point M~ and 
edge toward the mesh point (j — 1/2, n). Also note that, by using the above preliminaries, one can show 
that: (i) 

r = 0 if v = Q (3.17) 

(ii) 

lim r = 1 (3.18) 

V 2 —>1~ 

and (iii) 

T > T 0 (v 2 ) (v 2 < 1) (3.19) 

Recall that (i) (u$)" = (u§)™ if r = v = 0; and (ii) (u$)” = (u|)” if r = 1. As such, Eqs. (3.17) and 
(3.18) imply that, for each member in the subset, (i) (u$)" = (u|)” if v = 0; and (ii) (ft$)" = (u|)" in 
the limit of v 2 — » 1~. In other words, all members in the subset are ideal solvers in the domain v 2 < 1. 
Moreover, by using Eq. (3.12) and (3.19), one can also show that these ideal solvers are also stable in the 
same domain. Hereafter, each of these ideal solvers will be referred to as an ideal c-r* scheme. 

Corresponding to infinitely many choices of h, there are infinitely many different ideal c-t* schemes. In 
particular, because 

h(s) = y/s (0 < s < 1) (3.20) 

is a legitimate choice (see Eq. (3.14)), Eq. (3.16) implies that an ideal c-t* scheme can be defined using the 
relation 

r = \v\ (v 2 < 1) (3-21) 

Any ideal c-r* scheme described above meets all the requirements of an ideal solver defined in Sec. 2.2. 
However, because of stability problem, an ideal c-t* scheme and its multidimensional extensions may not 
be robust enough for some complicated real-world applications. To overcome this difficulty, we introduce a 
special c-r scheme which is defined by Eqs. (2.8), (3.9), and (3.10) with 

t = P M (P > 1; J' 2 < 1) (3.22) 
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where (3 > 1 is an adjustable parameter. For this scheme, Eq. (3.17) is also valid. Thus it becomes the 
nondissipative a scheme when v = 0 and, therefore, it will not become overly dissipative as v — > 0 (i.e., the 
scheme’s numerical dissipation is Courant number insensitive). Moreover, because (3 > 1, Eqs. (3.12) and 
(3.14) along with the comment (b) given following Eq. (3.14) imply that the scheme is stable if v 2 < 1 and 
it becomes more dissipative (i.e., more stable) as the value of (3 increases. However, also because (3 > 1 , 

lim t = (3 > 1 (3.23) 

bhi- 


As such, except for the case (3 = 1, the scheme does not become the c scheme in the limit of \is\ — > 1”, i.e., 
it does not meet all the requirements of an ideal solver. 

Note that, unless specified otherwise, in Sec. 4 we consider only the cases in which either Eq. (3.16) or 
Eq. (3.22) is assumed. 

4. Extensions of the c-t* scheme and related weighted averagings 


To proceed, let 


and 


(fix-)? 


def AX 

= T 


( «? - u(P~) \ 
y (1 + t) ax/4 J 


(fix+)? 


def AX 

“ T 


( u(P + )-v? \ 
y (1 + r ) ax/4 J 


(4.1) 


(4.2) 


Because |AP - | = |AP+| = (1 + t)ax/4 (see Fig. 3), it is easy to see that (u $ _)? and (u$+)" are two 
one-sided difference approximations of du/dx at the mesh point (j, n) with one being evaluated from the left 
and another from the right. Also, it follows immediately from Eqs. (3.8), (4.1) and (4.2) that 


(fix)? = \ [(fix-)? + (fix+)?] 


(4.3) 


Moreover, by (i) substituting Eqs. (2.8), (3.6) and (3.7) into Eqs. (4.1) and (4.2), and (ii) using Eqs. (3.10) 
and (3.17), one arrives at the conclusion that 


(fix-)? = (fi*+)? = (fix)? iy = 0) 


(4.4) 


when v = 0. 

With the above preliminaries, several extensions of the c-t* scheme will be constructed in the following 
subsections. 

4.1. Scheme w-1 

A comparison of Eqs. (4.1)-(4.3) with Eqs. (2.32)-(2.34) reveals that an obvious extension of the c-t* 
scheme can be obtained by replacing (u $ _)? and (u$+)? in Eqs. (2.35) and (2.39) with (fig-)? and (u £+ )?, 
respectively. In other words, the new extension is formed by Eq. (2.8) and 

(«x)? = (w-)?(fi,-)? + ( u; +)?(fix+)? (4.5) 

with 

M] = W±((fi*-)?, (fix+)?, «) (4.6) 

Because the scheme is the first extension of the c-t* scheme in which (it®)? is expressed as an weighted 
average of (ftg_)? and (u $+ )?, for simplicity, hereafter it will be referred to as Scheme w-1. It has been 
shown numerically that Scheme w-1 is stable if \v\ < 1 and a > 0. 

Note that, as a result of Eqs. (2.38), (4.4) and (4.6), one concludes that, for any given a > 0, (w_)? = 
(tu+)? = 1/2 if v = 0. In other words, for Scheme w-1, the “weighted” average on the right side of Eq. (4.5) 
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becomes a simple average if v = 0. According to an explanation given in the last paragraph of Sec. 2, this 
implies that Scheme w-1 will lose its capability to suppress wiggles or overshoots when v becomes small. 
For this reason, even though the Euler version of Scheme w-1 performs much better than that of the special 
scheme referred to in Sec. 2.3 in its ability to resolve shocks and contact discontinuities crisply in a wide 
range (from 1 to less than 0.001) of the global Courant number (i.e., the maximal value of local Courant 
numbers), it has a serious shortcoming, i.e., wiggles or overshoots can appear near a discontinuity in a 
generated solution when the local Courant number there becomes extremely small. In the following, it will 
be shown that this weakness can be overcome by simple modifications of Eq. (4.6). 

4.2. Scheme w-2 

A new scheme, referred to as Scheme w-2 is formed by Eqs. (2.8) and (4.5) with (w±)" being given 
by Eq. (2.39). In other words, although (it $ )” is still constructed as an weighted average of (tt $ _)™ and 
(v,x+)'j, the associated weight factors (■ w±) 1 j are evaluated using («$-)" and (% + )". Because the last two 

parameters, respectively, are identical to the special cases of (u s _)" and (u$+)” with r = f 1 (see Eqs. (2.12), 
(2.32), (2.33), (3.4), (3.5), (4.1), and (4.2)), their values do not vary with v. As such, [w±) r - ^ 1/2 and 
therefore the weighted average on the right side of Eq. (4.5) will not turn into a simple average when v = 0. 
In other words, Scheme w-2 is still capable of annihilating the numerical wiggles near a discontinuity even if 
v becomes small. It has been shown numerically that Scheme w-2 again is stable if \v\ < 1 and a > 0. 

Note that a possible drawback of Scheme w-2 is that the relation |(u s _)"| < |(u $+ )"| (|(it s _)"| > 
|(tix+)"|) does not automatically follow from |(tt $ _)”| < |(u $+ )”| (|(u s _)"| > |(u $+ )"|) and vice versa. 
As a result, at some local mesh points, it may happen that, of («i_)" and (%+)", the one with smaller 
absolute value may not be associated with a weight factor > 1/2. According to a discussion given in the 
last paragraph of Sec. 2, this implies that there is no guarantee that, at all localities, the weighted-averaging 
induced numerical dissipation will be available to suppress wiggles or overshoots. Despite this possible failing, 
fortunately it has been demonstrated numerically that, not only are they capable of suppressing wiggles or 
overshoots robustly, Scheme w-2 and its Euler extensions are also highly accurate. 

In the following, schemes that overcome the weakness of Scheme w-1 and also avoid the theoretically 
possible failing associated with Scheme w-2 will be constructed using new weighted-averaging formulae more 
advanced than that given in Eq. (2.38). 

4.3. New weighted-averaging techniques 

To pave the way, first we shall discuss a limitation of Eq. (2.38) as a generator of weight factors. 

Let x± ^ 0. Then, for a given a > 0, obviously W- — > 1/2 and W + — > 1/2 as \x + /x- \ — * 1. As such, 
when \x+/x- \ is very close to 1, then both W- and W+ will be very close to 1/2 unless a 1. As a result, 
in case that (i) (ux±)*j yf 0, (ii) |(u $+ )”/(u $ _)"| is very close to 1; and (iii) Eqs. (4.6) is assumed, then the 
only way to prevent the weighted average that appears on the right side of Eq. (4.5) from becoming almost 
a simple average is to increase the value of a used. However, this approach may be impracticable because 
numerical evaluation of a quantity such as x a for any real number x generally is hampered by round-off 
errors and thus becomes highly inaccurate if the value of a becomes too large, say 100. It is the purpose 
of this subsection to introduce new weighted-averaging techniques that do not have the limitaton discussed 
above. 

For motivation, note that Eqs. (4.5) and (4.6) can be expressed as 



('Ux)j = W\X\ + W 2 X 2 

(4.7) 

and 

w\ = — — — and w 2 = — — — (si + S 2 > 0) 

Si + S 2 Si + S 2 

(4.8) 

respectively if 

Xi d = (Ux-)j and x 2 d = ( Us+)j 

(4.9) 


wi = f ( w _)" and w 2 = f (tu+)" 

(4.10) 
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Si = |(*W)"r and s 2 = !(«*-)”!“ (a > 0) (4.11) 

Equation. (4.7) represents an weighted average of only two values X\ and X 2 ■ However, for the sake of 
generality, weighted averages of two or more values will be considered in the following development. 

To proceed, let (i) N be an integer > 2, (ii) se, £ = 1, 2, . . . , TV, be given positive numbers, and (iii) 


def S£ 


e = i,2,...,N 


(4.12) 


where 

»= (l>) >o 


(4.13) 


(Note: to streamline the following development, here we assume that se > 0, £ = 1, 2, . . . N, instead of se > 0, 
£ = 1, 2, . . . , N, as could be inferred from Eq. (4.11). However, without causing any practical harm, one can 
add a very small positive number, such as 10“ , to each member of a set of nonnegative numbers and turn 
all of them into positive numbers). It follows from Eqs. (4.12) and (4.13) that 


N 

we = 1, and 1 > we > 0, £ = 1,2, . . . , N 

e=i 


(4.14) 


As such, 


N 


w = 


we xi 


(4.15) 


1=1 

is an “interpolated” weighted average of the real numbers xe- Note that, unless specified otherwise, hereafter 
i = 1, 2, . . . , N is assumed. 


Let 


Then 





(4.16) 

(4.17) 


Also, with the aid of Eq. (4.14), Eq. (4.16) implies that 

N 

J2 s c = ° ( 4 - 18 ) 

1=1 


Note that W becomes the simple average of xe, £ = 1, 2, . . . , N, if all Se = 0. Thus the set {(5i, A 2 , . . . , Sn} 
provides a measure of how far the weighted average is deviated from the simple average. In the following, a 
simple way to adjust this deviation will be introduced. 

Let 

'Wn = f min {M and Vax = f max{<^} (4.19) 

Then Eq. (4.17) and the fact that 1 > we > 0 for all £ imply that 

1 > — + <5max and — + <5 m j n > 0 (4.20) 

Let some Se ^ 0 (i.e., the case with all we = 1/iV is excluded). Then Eq. (4.18) implies that imax > 0 > t> m ; n . 
The last inequality and Eq. (4.20) can be combined to yield 

1 - If > <5max > 0 > <5 min > -4 (some Se ^ 0) (4.21) 
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Note that an immediate result of Eq. (4.21) is 


def 

v max = mm 


Umax 




N5„ 


> 1 


(some 5( ^ 0) (4.22) 

Given any adjustable real parameter a > 0, let 

5’ e d =a6 e (4.23) 


Then Eq. (4.18) implies that 


In turn Eq. (4.24) and 


imply that 


N 


E^ = ° 


e=i 


/ def 1 

1 = N + ° e 


N 


E w ^ = i 


t = i 


(4.24) 


(4.25) 


(4.26) 


As such, w' e , £ = 1,2 , . . . ,N, form a new set of weight factors. From Eqs. (4.23) and (4.25) one also concludes 
that the disparity of the weight factors (i.e., the deviation of the values of the weight factors from 1/N) will 
be amplified (reduced) if a > 1 (a < 1). 

The condition that 

1 > v/f > 0 (4.27) 

will be imposed in the current development. With the aid of Eqs. (4.19), (4.23), and (4.25), and the original 
assumption that a > 0, one concludes that Eq. (4.27) is true if and only if 

1 > — + a <5max and + a ^'min — ® (4.28) 

For the special case that in which all Si = 0, one has <5 m ; n = <5 m ax = 0 and thus Eq. (4.28) is true always. 

On the other hand, for the case that some Si 0, Eqs. (4.21) and (4.22) imply that Eq. (4.28) is equivalent 

to 

Umax > a (some Se ^ 0) (4.29) 


Let some Si ^ 0. Then, according to Eq. (4.22), er ma x > 1- Moreover, Umax increases as |<5max| and 
(5 in inl decrease. In fact, Umax — > +oo as <5max — * 0 + and 5 m j n — > 0 _ . Thus the range of the values of a 
allowed becomes larger when |dmax| and <5 m ; n become smaller. Note that, when W defined in Eq. (4.15) 
almost becomes a simple average (i.e., when |£ m ax| 1 and |<5 m j n | <C 1), the disparity of the weight factors 
must be sharply amplified such that the weight average 


W' 


def 


N 

xt 


t = i 


(4.30) 


will deviate substantially from the simple average. In this case, the larger range of the values of a allowed 
meets the need to use a larger value of a. In practice, the value of a used can be generated using a preset 
formula as long as the generated value is less than or equal to Umax- For the case that the value generated 
using the preset formula is larger than cr m ax, cr = Umax is assumed. 
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As an example, consider the N — 2 case in which xg and sg, £ = 1,2, are defined by Eqs. (4.9) and 
(4.11). It was explained earlier that, for this case, w\ — > 1/2 and w 2 — > 1/2 as v — » 0. In other words, 
the weighted average wix± + W2X2 almost becomes a simple average when \v\ <C 1. To prevent this from 
happening, the weight factors w± and W2, respectively, are replaced by the new weighted factors w[ and w 2 
generated assuming 


a = min 


Umax, 



(cr 0 > 0) 


(4.31) 


where cr 0 > 0 is a preset parameter in the order of 1. According to Eqs. (4.16), (4.19), and (4.22), the fact 
that w\ — > 1/2 and W2 — ► 1/2 as v — > 0 implies that cr m ax — ► 00 as v —> 0. As such Eq. (4.31) implies that 
cr 1 when \v\ -C 1. 

Note that, for any N = 2 case, one of 8\ and 82 is <5 m ax while another is <5 m j n . As a result, Eqs. (4.18), 
(4.20), and (4.22) imply that 


0 < <5max — — ^min ^ 1/2 (some 8g ^ 0) (4.32) 

and 

Umax = — x (some 8 g 0) (4.33) 

& <>max 

Also for any case with N = 2, 5 m ax > 0 and a = crmax, Eqs. (4.23), (4.25), (4.32), and (4.33) imply that: 
(i) w[ = 1 and w' 2 = 0 if <5i = <5max, and (ii) w' 2 = 1 and w[ = 0 if 82 = ^max- 

This completes the description of a new approach by which the weight factors tu/, i = 1,2,..., TV, are 
generated from the given weight factors wg, £ = 1, 2, . . . , N. In the following, Another approach will be 
described. 

To proceed, the indices of sg, £ = 1, 2, . . . , N, will be reshuffled such that 


sn > sjv _ 1 > . . . > Si > 0 (4.34) 

As such, Eqs. (4.12) and (4.13) imply that Eq. (4.14) can be replaced by a set of stronger conditions, i.e., 


N 


wg = 1 and 1 > wn > Wn - 1 > . . . > Wi > 0 




Next let 


def S£+ 1 

m = 1 , 1 = 1 ,.. • , N-l 


sg 


(4.35) 


(4.36) 


Then (i) ijg > 0, £ = 1, . . . , N — 1, and (ii) 


sg+i — 


n (!+%') 


'=1 


si, 


£=1,...,N-1 


Given any adjustable real parameter a > 0, let (i) §1 = si and 


(4.37) 


and (ii) 


sg + 1 — 


" e 

11(1 + "*) 


U'=l 


si, 


£= 


~ def Sg 
Wg = — , 

S 


£=1,2,.. .,7V 


(4.38) 


(4.39) 
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where 


s = (i>) >° 

Because cr > 0 and rje > 0, £ = 1, . . . , IV— 1, Eq. (4.38) implies that 


(4.40) 


SjV -2 SjV— 1 >...>§ 1 > 0 


(4.41) 


Also, as a result of Eqs. (4.39)-(4.41), one has 

N 

we = 1 and 1 > wn > u>n-i > . . . > U)i > 0 (4.42) 

e=i 

As such, we, £ = 1,2 ,N, form a new set of weight factors and 

N 

W = f ^ we x e (4.43) 

e=i 


is an “interpolated” weighted average of the real numbers xe- Note that, for the special case that sn = sjv-i = 
. . . = si > 0, it is easy to see that (i) we = we = 1/N, £ = 1, 2, . . . , N, and (ii) rje = 0, £ = 1, . . . , N — l. 

Let l\ and £2 be any pair of integers with 1 < £\ < £2 < N. Then Eqs. (4.12) and (4.37)-(4.39) imply 


that 


and 


nV+’fcO 

Wo, - LA 

£'=h 

(4.44) 

e 2 -i 

^ = ri( i +^') 

Wo, 

£'=ix 

(4.45) 


Because cr > 0 and rje > 0, £ = 1, . . 
we 2 /u>e 1 = 1 if i)e = 0 for all £ with 
£1 < £ < {£2 — 1), one has 


. , N—l, a comparison of Eqs. (4.44) 
£\ < £ < (£2 — 1). However, in case 

f > ^ if cr > 1 

w ‘i 


and 

that 


(4.45) reveals that 
rje 7^ 0 for at least 


we 2 /w tl = 

one £ with 


we 2 

Wti 


Wl, 

we. 


if cr < 1 


w e 2 

we. 


if cr = 1 


(4.46) 


^From the above discussions, one concludes that, except for the special case in which sn = sjv-i = . . . = 
Si, the disparity of uie is greater (less) than that of if cr > 1 (cr < 1). Note that the current approach 
for amplifying the weight factors has one advantage over the approach described earlier, i.e., in the current 
approach, there is no upper bound for the value of a one could use. Thus, in the current approach, Eq. (4.31) 
can be simplified as 

^ = n ( 4 - 47 ) 

M 

where cr Q > 0 again is a preset parameter in the order of 1. 

Note that, after sorting through the differences in the notations used in [49] and here, and using the 
fact that the value of the smaller of the parameters (s_)" and (s+)" defined in Eq. (3.23) of [50] is zero, it 
can be shown that the weighted-averaging technique introduced in Eqs. (3.23), (3.26), and (3.27) of [50] is 
equivalent to the special case N = 2 and cr Q = 1/2 of the second approach just described above. 
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4.4. Schemes w-3 and w-4 

Consider the N = 2 case in which xg and se, i = 1,2, are defined by Eqs. (4.9) and (4.11). Let (' w'_) f 
and (w '_ |_)", respectively, be the weight factors associated with (u s _)" and (m s+ )" generated using the first 
approach described in Sec. 4.3. Then, by definition, Scheme w-3 is formed by Eq. (2.8) and 

(«*)? = («/_)7(«a_)7 + K)?(u a+ )? (4.48) 

On the other hand, let (u)_)" and (w + )™, respectively, be the weight factors associated with (u s _)" 
and (u s+ )” generated using the second approach described in Sec. 4.3. Then, by definition, Scheme w-4 is 
formed by Eq. (2.8) and 

(«*)? = («)-)?(«*_)? + (tD + )?(u* + )? (4.49) 


5. Courant number and Mach number insensitive solvers 

We consider a dimensionless form of the 1-D unsteady Euler equations of a perfect gas. Let p, v, p , and 
7 be the mass density, velocity, static pressure, and constant specific heat ratio, respectively. Let 

«i = P, u 2 = pv, u 3 = p/( 7 - 1) + (1/2 )pv 2 (5.1) 

fi = u 2 (5.2) 

h = (7 - 1)«3 + (1/2) (3 - "i){u 2 ) 2 /u\ (5.3) 

and 


h = iu 2 u 3 /ui - (l/2)(7 - 1 )(k 2 ) 3 /(mi ) 2 
Then the Euler equations can be expressed as 


(5.4) 


du r 


dfn 


dt dx 

The integral form of Eq. (5.5) in space-time E 2 is 


= 0, m = 1, 2, 3 


(5.5) 



• ds = 0, 


in = 1,2,3 


(5.6) 


where h m = (/ m , u m ), m = 1,2, 3, are the space-time mass, momentum, and energy current density vectors, 
respectively. 

As a preliminary, let 

fm,k = f df m /duk, m,k = 1,2,3 (5.7) 

Let F be the Jocobian matrix formed by f m ,k , m,k = 1,2,3. Then Eqs. (5.2)-(5.4) imply that 


F = 


7-3 ( u 2 




(3-7) 


U2 

Ml 


M 2 M 3 M3 _ 3(7 - 1 ) / U 2 

(mi) 2 7 m 1 2 


0 \ 

7-1 


/— 

Ml 


J 


Let c be the sonic speed. Then 


c = 


7(7 - 1) 


M3 

Ml 


1 

2 



(5.8) 


(5.9) 
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Let G be the 3x3 matrix defined by 


G 


def 


1 

V 2 


U2 

Ml 

U2 

Ml 


M 2 

2 i/ 2 cmi 


m 1 

V2c 

M 2 Ml 

Tic Ti 

T2 


MiC 


M 2 

1 / 2(7 — 1 ) 2 -\/ 2 cmi 


Ml 
V2 c 

m 2 mo_ 

T2C i/2 

-ui*.' 

T2 


(5.10) 


MlC 


T2(7 - 1) / 


Then the inverse of G is given by 


G~ x 


( 


\ 


(7- l )( u 2 ) 2 
2 V / 2 c(mi ) 3 
(7- 1)(^2) 2 
2 V / 2 c(mi ) 3 



m 2 

T2(Mi) 2 


M2 

T 2 (Mi) 2 


(7 - 1 )m 2 1 — 7 \ 

C 2 M1 c 2 

(1 — 7)m 2 1 7 - 1 

-\/2c(mi) 2 \/2mi v^cmi 

(1 - 7)m 2 1 7 - 1 

i/2c(mi) 2 \/2mi y / 2 cui / 


(5.11) 


Moreover, for any numbers a±, o 2 , . . ., a n , let diag(a \ , o 2 , . . . , a n ) denote the diagonal matrix with ai, o 2 , 
. . ., a n being the diagonal elements on the first, second, . . ., and n-th rows, respectively. Then, by using 
Eqs. (5.8)-(5.11) and v = u 2 /u\, one has 


G 1 F G = diag(v, v — c, v + c) 


(5.12) 


For any (j,n) € and any (x,t) € SE u m (x,t), f m (x,t), and h m (x,t) are approximated by 
m^(x, t ; j, n), f^(x,t ; j,n), and h* m {x,t\ j,n), respectively. They will be defined shortly. Let 


u m( x , t ; j, n) = {u m )j + {u mx )j{x - Xj) + 


t n ) 


(5.13) 


where (M m )”, (u mx )J, and (M mt )" are constants in SE(j, n). 

By definition, for each m = 1,2,3 and each k = 1, 2, 3, f m and / m> fc are functions of Mi, m 2 and M3. Let 
(/ m )” and ( denote the values of f m and f m ,k , respectively, when the independent variables Mi, m 2 , 
and M3, respectively, assume the values of (mi)”, (m 2 )", and (M3)”. Let 


(fr, 


0? = y!(/ n 


fc=l 




and 


Because 


and 


UmtTj = £(/m, *)>**)" 

k = 1 


— ^ , fm,k 


du k 


dfrr 

dx ~^ jm ’ K dx 

fc= 1 


5/m _ , dMfc 

dt - 2 ^^^ 

fc= 1 


(5.14) 


(5.15) 


(5.16) 


(5.17) 
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(/„)“ and (/mt)" can be considered as the numerical analogues of the values of df m /dx and df m /dt at 
(. Xj,t n ), respectively. As a result, we assume that 


/mOM 5 d = (fm)j + ( fmx)](x ~ Xj) + ( fmt)j(t - t n ) 
Because h m = (f m , u m ), we also assume that 


t ; j, n ) = f (/^(x, t ; j, n), <,(x, t ; j, n)) 


(5.18) 


(5.19) 


Note that, by their definitions, (i) (/ m )" and (/ TO ,fc)”, m = 1,2,3, are functions of (u m )", ur = 1,2,3, (ii) 
m = 1,2,3, are functions of ( u m )j and (w mx )", m = 1,2,3, and (iii) are functions of (u m )J 

and ( u mt )j , m = 1,2,3. 

Moreover we assume that, for any (x, t) € SE(j, n), u m = u^Xjt ; j,n) and / m = /j^(x,t; j, n) satisfy 
Eq. (5.5), i.e., 


du* m (x,t-,j,n) df^(x,t\j,n) 
dt dx 


(5.20) 


According to Eqs. (5.13) and (5.18), Eq. (5.20) is equivalent to 


(■ Umt )7 


-{fmxYi 


(5.21) 


Because (/ TOX )" are functions of (u m )" and (u mx )", Eq. (5.21) implies that ( u m t )" are also functions of 
( u m )” and (u mx )”. From this result and the facts stated following Eq. (5.19), one concludes that the only 
independent discrete variables needed to be solved in the current marching scheme are (rt m )" and («„)". 
In the current development, the Euler counterpart of Eq. (2.11), i.e., 


/ 

Js{CE(j,n )) 


• ds = 0 


is assumed. For any (j,ri) £ fi, let 


— f 

V 11 mx ) j — 4 v 11 mx ) j 


and 


/ \n i / j? \ 7i | (^0 / p \n 

{Sm)j — {Umx)j + + 4 Aa . \fmt)j 


Then, with the aid of Eqs. (5.13)-(5.15), (5.18), (5.19), and (5.21), Eq. (5.22) implies that 


(u m Y = 


1 

2 L 


( u '| n - 1 /2 , / ^-1/2 , / \n-l/2 _ / 1/2' 


(5.22) 


(5.23) 

(5.24) 


(5.25) 


Eq. (5.25) forms the first component of each of the Euler schemes to be constructed here. As will be shown, 
the second component which evaluates (it mx )" is scheme dependent. 

5.1. The Euler c-r scheme 


To proceed, consider any ( j , n) 
matrices formed by u m , ( u m )", (it. 
v, the sonic speed c and the matrices F, G, and G _1 

n n pn 
u 3 ’ 3' 3 ’ 

variables it, 
let 


€ ft. Let ft, Uj, (u x )", (it x )”, and (u t )", respectively, denote the column 
TOX )", (u TOX )", ( Umt . )", m = 1,2,3. Moreover, recall that the fluid velocity 


are functions of u m , m = 1,2,3. Thus we will define 
G", and (G -1 )™, respectively, to be the values of u, c, F, G, and G -1 when the independent 
re , m = 1, 2, 3, respectively, assume the values of ( u m )!?, m = 1, 2, 3. Given the above definitions, 


= f ^ 


At 


AX 


(^)? = 


def (U - c)”At 


AX 


and („)» S 

J AX 


(5.26) 
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and 


(5.27) 


(*w)? = (Kl + |c?l)f- 

j j j 

Because 

{vmaxYj = maxiK^)"!, K^)"!, |(i / 3)j |} (5.28) 

( Vmax )™ can be interpreted as the Courant number at the mesh point ( j,n ). 

Next, for each m = 1,2,3, let the points P+ and P~, and the parameter (r m )” shown in Fig. 5 be 
defined in the exact same manner by which the points P + and P ~ , and the parameter r were defined (see 
Fig. 3). Moreover, let (i) [b] m denote the rrath component of any column matrix b, (ii) 


[(G- 1 )"u] m (P+) d ^ f (G _1 )"u"“jy2 + (At/2) -[1 (G- 1 )]^)^ 


-m- 1/2 


n-1/2] 


\n— 1/2 


and (iii) 


(5.29) 


[(G- 1 )]u] m (p-) d ^ [(G- l )^: 1 1 / 2 2 ] m + (At/2) [(G- l )”(u,);r 1 % 2 ] m + [i-(r m )^ 

(5.30) 

Note that, in the current development of the procedures for evaluating (u s )", (G -1 )!/ is treated as a fixed 
constant square matrix for any given fixed (j, n) € Q while u is treated as a variable column matrix (this 
practice, in spirit, is similar to the definition of u*(x,t;j,n ) given in Eq. (2.3) where w”, ( u x )", and (iq)” 
are treated as constants while x and t are treated as variables). Thus 


8 [(G _1 )"5 


dt 


(G Gdi 


and 


d [(G _1 )"S 
dx 


(G- 1 ) 


, du 


— 1 \n 


(5.31) 


It follows that 


(G _1 )j (^)”±i/2 


and 


(G 1 )"(ux)” ± i 1 /2 > respectively, are the numerical analogues of 


d [(G 1 )j"u] m /dt and d [(G 1 )"u] m /9x at the mesh point (j ± 1/2, n — 1/2). With the aid of this interpre- 
tation, Eqs. (5.23), (5.29), and (5.30) imply that [(G -1 )"u] m (P+) is a first-order Taylor’s approximation of 
[(G- 1 ) j"u\ at point P+ evaluated using the marching variables at the mesh point ( j + 1/2 , n — 1/2) while 
[(G- 1 ) j"u\ m (P m ) JS a first-order Taylor’s approximation of [(G 1 )"u] m at point P m evaluated using the 
marching variables at the mesh point (j — 1/2, n — 1/2). As such, Eqs. (5.29) and (5.30) can be considered as 
the Euler versions of Eqs. (3.4) and (3.5), respectively. In the following, we will construct the Euler versions 
of Eqs. (3.6)-(3.9). 

By using Eqs. (5.14), (5.21), and (5.23), one has 


(At/2)(G~ 1 )](u t )]~lG = -(2At/Ax)(G- 1 )]F^GG](G~ 1 )](u^~l / ^ 


Let F™±i /2 ~ P/ ■ Then Eqs. (5.12) and (5.26) imply that 


At ^- h „ P n-i/2^ ^ - ( G-i)"F"G" = diag((v^,(v^,(v^) 


( rt— i\n 771/1— 1/ z 


AX J J 1 J AX 

Combining Eqs. (5.32) and (5.33), one has 

t —1 \ nf^t \ n — 1/2 


(At/2) (G- 1 )?^);-^ « -2(v m )3 I (G~ L )™ (uxYj±ljl 


(5.32) 


(5.33) 


(5.34) 


Substituting Eq. (5.34) into Eqs. (5.29) and (5.30), one arrives at the current Euler versions of Eqs. (3.6) 
and (3.7), i.e., 


[(G- l )-u\ m (P+) « [(G- 1 )”^ 1 / 2 1 - [2{v m yy + 1 - (r m )]] [(G-^^T-ljl 


- 1/2 


(5.35) 
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and 


[(G~ 1 )]u] m (P-) « \{G- 1 )^Zl / / l] m [2(^m)" - 1 + (r m )”] [(G-^iu^H 

Next we introduce the Euler version of Eq. (3.8), i.e., 


(5.36) 




def AX [(G~ 1 )?u] m (P+) - [(G-^u] m (P~) 


[l + (T m )j] AX/2 


((r ro )^-l), m= 1,2,3 (5.37) 


where [(G 1 )”u 1(P+)and [(G- 1 ) "u] (P m ), respectively, are to be evaluated using Eqs. (5.29) and (5.30) 

instead of the approximated formulae Eqs. (5.35) and (5.36). With the aid of Fig. 5, Eq. (5.37) implies that 
[(G- 1 ) j’u] - (j, n) is a central-difference approximation of d [(G 1 )"u] m /dx at the mesh point (j, n). 

To proceed, note that 


Thus 


d[(G^u] = . M 

dx 1 >j dx 


du d[(G~^u] 

dx J <9x 


Let [(G 1 )”m] 2 (j, n) be the column matrix formed by [(G *) 
Eq. (5.37)). Then, with the aid of Eq. (5.39), one concludes that 


~ L jU\ _ 
J -i mx 


(5.38) 

(5.39) 

( j,n ), to = 1,2,3 (which are defined in 


(4)? = ^{[(G-^.y.n)} 

represents a central difference approximation of du/dx at the mesh point (j, n). Thus 

(**)? = (4)” 


(5.40) 


(5.41) 


is an Euler version of Eq. (3.9). The Euler c-r scheme is defined by Eqs. (5.25), (5.29), (5.30), (5.37), (5.40) 
and (5.41). It has been shown by numerical experiments that stability of this scheme generally requires that 

{"max)] < 1; and (r m )" > To (((r'm)") 2 ) ((j,n)ef2), to =1,2, 3 (5.42) 

Obviously, the form of Eq. (5.42) is very similar to that of Eq. (3.12). 

Consider the special case in which (n)" = (t 2 )" = ( 73 )" = r”. For this case, (i) points P^, P^ and P 3 + 
coincide (the resulting common point will be denoted by P + ); and (ii) points Pj - , P 2 _ , and P 3 _ also coincide 
(the resulring common point will be denoted by P~). As such it can be shown easily that Eqs. (5.29) and 
(5.30) reduce to 

u(P + ) = [rt+ (Ai/2)t? t - (1 - tJ)Ux\\ 

and 


respectively. Also Eqs. (5.37) and (5.40) can be used to conclude that 


(fix)" = 


def u(P + ) - u(P ) 


2(1 + t") 


n- 1/2 
J+l/2 

(5.43) 

n- 1/2 

J — 1/2 

(5.44) 

-1) 

(5.45) 


Note that: (i) with the aid of Eq. (5.28) and the fact that r Q (s ) is strictly monotonically increasing in the 
interval 0 < s < 1, one concludes that Eq. (5.42) reduces to 


(y max) 1 / < 1; and t” > To (((rwx)") 2 ) (Cb n ) e 


(5.46) 
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for this special case; and (ii) Eqs. (3.4), (3.5), and (3.8), respectively, are very similar to their Euler versions 
Eqs. (5.43)-(5.45). 

The above simplified scheme has the advantage in its structural simplicity. However, it suffers from the 
loss of the capability to adjust the value of individual (r m )". This is a serious disadvantage because the 
theory presented in Sec. 3 along with the fact that Eqs. (3.6) and (3.7), respectively, are very much similar 
to Eqs. (5.35) and (5.36) in algebraic form strongly suggests that the Euler c-t scheme will perform better 
if, for each m, the value of (r m )" can be adjusted to match the value of t 0 (((z'm)") 2 )- This is particularly 
true if, for some (j. n) € fi, there is a large disparity among the values of |(^ m )"|, m = 1,2, 3, i.e., 

(w)? » min{|(zq)"|, |(^|} (5.47) 

As an example, consider an Euler flow solution with an extremely small Mach number everywhere (i.e., 
|u"| <C |c"|, ( j,n ) € fi), For such a case (see Eqs. (5.26)-(5.28)), 

(C i,n)en) (5.48) 

As such, Eqs. (3.13) and (5.46) can be used to show that 

T o (((Zff)") 2 ) « Tj (5.49) 

Because the c-t scheme tends to become more dissipative as the value of r increases from t 0 (v 2 ) (see comment 

(b) given following Eq. (3.14)), the fact that Eqs. (5.35) and (5.36), respectively, are very much similar to 
Eqs. (3.6) and (3.7) in algebraic form strongly suggests that the component [(G , ^ 1 )"u] 1 will become highly 
dissipative under the condition Eq. (5.49). As such , one expects that a solution to the simplfied Euler c-t 
scheme would become highly dissipative at a region with a small Mach number. 

5 . 2 . The Euler c-t* schemes 

Let h m (s ) (0 < s < 1), m = 1, 2, 3, be strictly monotonically increasing smooth functions which satisfy 
the current version of Eq. (3.15), i.e., 

h m (0) = 0; lim h m {s) = 1; and h m (s) > t 0 (s) if 0 < s < 1, to =1,2, 3 (5.50) 

S — > 1 — 

Then, based on the similarity in form that existed between Eqs. (3.6)-(3.8) and Eqs. (5.35)-(5.37), an ideal 
Euler c-t* scheme can be formed as a special Euler c-t scheme in which 

(TV*)? = h m ((KO?) 2 ) (KO?| < 1), TO = 1,2, 3 (5.51) 

Obviously, Eq. (5.51) is the Euler version of Eq. (3.16). 

For the same reason that justifies the use of the relation Eq. (3.22), a special Euler c-t* scheme can be 
defined as a special Euler c-t scheme in which 

(r m )” = (/U?|(<)?| ((/U"> 1;|(0?|<1), rn = 1,2,3 (5.52) 

Here (/3 m )” > 1, to = 1, 2, 3, are adjustable parameters which may vary from one mesh point to another. 

Corresponding to the simplified Euler c-t scheme defined by Eqs. (5.25), (5.41), and (5.43)-(5.45), one 
can construct an Euler c-t* scheme in which 

rj* = h ((K*a,)”) 2 ) (|(*W)?| < 1) (5.53) 

where h(s ) (0 < s < 1) is any strictly monotonically increasing smooth function which satisfies Eq. (3.15). 
Obviously, we can also construct another scheme in which 

Tj* = 1 % I (tw)? I (Pj > 1; I ( w)? I < 1) (5.54) 
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Here /?” > 1 is an adjustable parameter which may vary from one mesh point to another. 

5.3. Other Euler extensions 

Let 


and 


Ffr -l r . , • , def AX - [(G~^u] m (P-) 

K G 4 [l + (r m )"] ax/4 

def AX [(G- l )W](P+)-[(G-^} m 




= 4 


[1 + (r m )!?] ax/4 


(5.55) 


(5.56) 


By definition, [(G 1 )"u] m __ (j, n) and [(G 1 )"u] m _ (j, n) are two one-sided difference approximations 
of d [(G~ 1 )j'u] m /dx at the mesh point (j, n) with one being evaluated from the left and another from 
the right. Let (t to )" be defined using either Eq. (5.51) or Eq. (5.52). Then, for each m, one can define 
[(G -1 )"u]™_ (j,n) to be an weighted average of and [(G~ 1 )]u] m _ + (j,n) using any 

weighted-averaging technique described in Sec. 4. As an example, an weighted average constructed using the 
second approach described in sec. 4.3 with N = 2 is given by 


def 


[( G “ 1 )" S ]L c?» = 


[1 + CTmVm- }j [(G ^u] ( j , n) + [1 + <J m Vm+]] [(G (j, n) 


Here 


and 


(Vm±)- = 



[2 H- cr m (? 7 m _ + 

[( G_1 )"«L*± (•?» 

Vm+)} ? 

min | 

[( G “^L*+C?>) 

, [(G -1 )"u] _ (j, n) 

L v 7 J J mx— / 

} 


- 1 > 0 


( \n def cr 0 

[ m)j - IW? 


(5.57) 


(5.58) 


(5.59) 


with cr 0 > 0 being a preset number in the order of 1. Note that: (i) for each (j, n) £ Cl, the value of the 
smaller of (rj rn +)j‘ and ( ij m _)" is zero; and (ii) to avoid dividing by zero, in practice, a small positive number 
such that 10 -60 should be added to each of the denominators that appear in Eqs. (5.58) and (5.59). 

Let [(G _1 )”i7] (j, n) be the column matrix formed by [(G -1 )”m] _ (j, n), m = 1, 2, 3 (which are defined 

in Eq. (5.57)). Then, with the aid of Eq. (5.39), one concludes that 


(■^i )j = f G” { [(G _1 )"u]™ (/, n) j 


(5.60) 


represents an weighted-averaging approximation of du/dx at the mesh point (j, n). Thus an Euler weighted- 
averaging c-r* scheme can be formed by Eq. (5.25) and 


(« = m 


(5.61) 


For the simplified case in which (n)!? = (r 2 )^ = (r 3 )™ = t”, Eqs. (5.55) and (5.56) are replace by 


, d., »7-“( p -> 

“ i + T n 


(5.62) 


and 


(**+)? = 


def U{P + ) - U7 


j 


1 + TV 


(5.63) 
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respectively. Let (i) r” be defined using either Eq. (5.53) or Eq. (5.54); and (ii) (u m ®±)" denotes the mth 
component of (u$±)"- Then, for each to, an weighted average of (u m j_)" and (u mS+ )" constructed using 
the second approach described in Sec. 4.3 is given by 


(v w -1" = 
v U'mx ) j 


def [1 + C u m *+)? + [1 + <T?(Vm+) 1 }] (Umx-Yj 


[2 + a(r) m - + Vm+)\j 


(5.64) 


Here 


and 


(Vm±)l = 


\(Umx±)j\ 

min{|(u m$+ )"|, |(u m£ _)"|} 


- 1 > 0 


(5.65) 


„ def 

CT i = 


q-p 

(^maa;)j 


(5.66) 


with a 0 being a preset number in the order of 1. Obviously a simplified Euler weighted-averaging c-r* scheme 
can again be formed using Eqs. (5.25) and (5.61) if (u™) 7 ) denotes the column matrix formed by (u^ s )!?, 
to = 1,2,3. 

Accuracy of two special Euler c-r* schemes will be evaluated in Sec. 6. The first is a simplified scheme 
in which we assume that 


|(Onax)?| < 1; r" = /3|l/„ 


and 


{ymaxYj 


(j, n) G Cl 


(5.67) 


Note that: (i) /3 > 1 and cr 0 > 0 are preset numbers in the order of 1; and (ii) the mesh point dependent 
parameter /3™ which appears in Eq. (5.54) is replaced by the mesh point independent parameter (3 here. 
Moreover, for a reason given immediately following Eq. (5.49), one would expect that a solution to such a 
simplified scheme becomes highly dissipative at a region with a small Mach number, i.e., the scheme is Mach 
number sensitive. Therefore, the first scheme will be referred to as the “Mach number sensitive scheme” in 
Sec. 6. Nevertheless this scheme generally is still Courant number insensitive. In fact the so called “new ” 
Courant number insensitive solutions presented in Figs. 4- 7 of [50] and Figs. 9 and 10 of [52] are generated 
using the current simplified scheme with (3 = 1.0 and cr 0 = 0.5. 

In the second scheme to be evaluated in Sec. 6, we assume that 


I iymax) 1 ^ I < 1; (r m )” = /Sml^m)?!; and 


/ _ \n def 


(( j,n)ea ), TO = 1,2, 3 (5.68) 


Note that: (i) (3 m > 1, to = 1,2,3, and cr 0 > 0 are preset numbers in the order of 1; and (ii) the mesh 
point dependent parameter (/3 m )" which appears in Eq. (5.52) is replaced by the mesh point independent 
parameter (3 m here. As will be shown in Sec. 6, not only is it Courant number insensitive, the second scheme 
is also Mach number insensitive. As such it will be referred to as the Mach number insensitive scheme in 
Sec. 6. 

6. Numerical results 

In this section, accuracy of the Mach number sensitive and insensitive schemes defined at the end 
of Sec. 5 will be evaluated by comparing their numerical solutions with the known analytical solutions of 
several shock tube problems. Without exception, (i) the spatial computational domaim which is defined by 
—0.505 < x < 0.505, is divided into 101 uniform intervals with ax = 0.01; and (ii) the specific heat ratio 
7 = 1.4. 

At any time t > 0, the exact solutions of a set of shock tube problems considered here are given by 


(p,P,v) 


(0.125,1.0, V0l4£) ifcc<v / (n4^ 
(10.0,1.0, V0T4£) if x > \/0T4£f 


(6.1) 
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Here, for each problem in the set, £ is a defining constant parameter. As such, for each problem, (i) p and 
v do not vary with x and t; and (ii) there is a contact (density) discontinuity which moves with the velocity 
= v. Obviously the discontinuity occurs at x = 0 when t = 0. Moreover, because 7 = 1.4, Eq. (6.1) implies 


that 



VTh2 = 3.3466 if 2; < y/OAA^t 


V0l4 = 0.3742 if 2; > V(hl4£f 


( 6 . 2 ) 


and 


. [ V0.0125£ = 0.1118£ if 2; < V0?\At;t 

M = - = l (6.3) 

C U if 22 > 


First consider the case with £ = 0.01. According to Eqs. (6.1) and (6.2), the flow field is characterized by 
extremely small values of M (0.001118 < M < 0.01) and a contact discontinuity moving with an extremely 
small velocity (= 0.00374). This problem, designated as Problem no. 1, is solved by the Mach number 
sensitive scheme assuming At = 0.0024, a Q = 2.0, and (3 = 1.0. It is also solved by the Mach number 
insensitive scheme assuming At = 0.0024 and a 0 = (3\ = fa = @3 = 2.0. Based on Eqs. (6.1) and (6.2), and 
the given values of £, ax, and At, it is estimated that the maximum local Courant number ecountered in 
each simulation is 0.804. At t = 60.0 = 25000At (i.e., y/0.14 £t == 0.2245), numerical values of pressure and 
velocity obtained from both simulations match the constant exact solution values to at least seven significant 
digits. As such no graphical comparisons of these numerical variables with their exact solution values are 
given. On the other hand, numerical values of density and Mach number at t = 60.0 are compared with 
the exact solution values in Fig. 6. It is seen that, for this problem characterized by extremely small values 
of M, the contact discontinuity is resolved much more crisply by the numerical values generated using the 
Mach number insensitive scheme than those generated using the Mach number sensitive scheme. Note that 
this difference in the schemes’ capability to resolve the contact discontinuity could become more pronounced 
if the chosen value of (3 is raised from 1.0 to 2.0, i.e., the value shared by (3\, @ 2 , and (3^. However, for the 
current case with very small values of (iq)”, reducing each value of u 0 and f3i from 2.0 to 1.0 will result in 
computational instabilty — obviously, at some ( j,n ), the assigned values of (a 1 )" and (n)" become too small 
to sustain stabilty. 

Next consider the case with £ = 10.0. For this case, the flow field is characterized by relatively large 
values of M (1.118 < M < 10.0) and a contact discontinuity moving with a relatively large velocity (= 3.74). 
This problem, designated as Problem no. 2, is solved by the Mach number sensitive scheme assuming 
At = 0.0012, o 0 = 2.0, and (3 = 1.0. It is also solved by the Mach number insensitive scheme assuming 
At = 0.0012 and o 0 = (3i = P 2 = (3s = 2.0. It is estimated that the maximum local Courant number 
ecountered in each simulation is 0.851. At t = 0.06 = 50At (i.e., y/0 T4£f = 0.2245), numerical values of 
pressure and velocity obtained from both simulations again match the constant exact solution values to at 
least seven significant digits. Moreover, as shown in Fig. 7, for this problem with relatively large values of 
M, the contact discontinuity is resolved crisply by both schemes. A comparion of these results and those 
shown in Fig. 6 reveals that the Mach number sensitive scheme is indeed Mach number sensitive while the 
Mach number insensitive scheme is indeed Mach number insensitive. 

The last shock tube problem to be considered is Sod’s problem [61]. For this problem, (i) at t = 0, 


(p,P,v) 


(1.0, 1.0, 0.0) if 2: < 0 
(0.125,0.1,0.0) if 2: > 0 


(6.4) 


and (ii) 0.0 < M < 0.929 for all x and t > 0. This problem, designated as Problem no. 3, is solved by the 
Mach number sensitive scheme assuming At = 4 x 10~ 6 and o 0 = (3 = 1.0. It is also solved by the Mach 
number insensitive scheme assuming At = 4 x 10~ 6 and o Q = (3\ = P 2 — /?3 = 1.0. The estimated maximal 
local Courant number encounted in each simulation is 0.00088. At t = 0.2 = 50000At, the computed solutions 
are compared with the exact solution in Fig. 8. In spite of the extremely small maximal local Courant number 
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encountered, it is seen that the numerical results generated by both schemes match very well with the exact 
solution. 

Problem no. 3 is also solved by (i) the Mach number sensitive scheme assuming Af = 4xl0 -3 and cr 0 = 
(3 = 1.0; and (ii) the Mach number insensitive scheme assuming At = 4 x 10~ 3 and a 0 = Pi = fa = /?3 = 1.0. 
The estimated maximal local Courant number encountered for each simulation is 0.88. At t = 0.2 = 50a1, 
the computed solutions are compared with the exact solution in Fig. 9. It is seen that the numerical solutions 
shown in Fig. 8 do not deteriorated much from the current results generated with a much larger maximal 
local Courant number. As such, both schemes are indeed Courant number insensitive. 

7. Conclusions and discussions 

Generally speaking, a stable numerical marching for a non-linear problem requires the presence of 
a sufficient amount of numerical dissipation. However, accuracy of the numerical results, especially for 
an unsteady problem, will suffer if too much numerical dissipation is present. As such, a careful control 
of numerical dissipation is a must for an accurate and stable non-linear unsteady numerical simulation. 
However, a proper control of numerical dissipation is a very difficult task. Although one can increase the 
numerical dissipation rather easily, it is much harder to reduce it when accuracy consideration requires it. 

The CE/SE method is developed from a set of non-dissipative solvers. As such each CE/SE solver is 
an extension of a core non-dissipative scheme. It is this unique feature that make it much easier to reduce 
numerical dissipation in a CE/SE simulation. It is also the key reason behind the successful construction of 
the Courant number and Mach number insensitive Euler solvers described in this paper. 

For the 2D and 3D unsteady Euler equations, there are two and three associated Jacobian matrices, 
respectively. In general, it is impossible to diagonalize these associated Jacobian matrices simultaneously 
using the same diagonalization matrix. Thus extension of the current work to a space of higher dimension 
is by no means trivial. 
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Figure 1 . — A surface element on the boundary 
S(V) of an arbitrary space-time volume V. 
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Figure 2. — The SEs and CEs. 
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